Modernization and Secularization

Does modernization lead to secularization in societies?

Let’s unpack some definitions:

Modernization refers, broadly, to economic development and institutions with some relation to economic development. Per capita GDP (PPP) is an obvious direct measurement of this, but it is also incomplete: a country may have a high average due to the presence of extractable natural resources that often lead to extremely skewed wealth distributions and, as such, averages that aren’t necessarily reflective of overall development. Life expectancy, safety standards, etc. also reflect modernization.

Secularization refers to a general shift away from religious practice and religiously-based values. This can be a bit nebulous as different faiths may emphasize different values, and is particularly problematic when comparing countries with different historical faith supersets (e.g. Islamic vs. Christian vs. Buddhist). For this analysis, I’ll include only countries in Western and Eastern Europe, North and Latin America, and Australia/New Zealand. While it would be beneficial to include sub-Saharan African countries, insufficient data exists in a number of public opinion-related questions, and as such, inference on them would rely heavily on imputation.

The basic idea of the modernization hypothesis is that as countries become wealthier, their citizens become more secular/less religious. The nature of this mechanism (by mechanism, we could essentially think of why it is that a maps to b) isn’t exactly universally agreed upon. For example, an economic theory might say that religion is a form of social insurance (churches often provide social services, and faith may in and of itself also act as an intangible insurance–a future payout for enduring hardships of today) and that when the state and industries within the state provide this additional insurance through better public and private goods, demand for the good diminishes. For our purposes, I’ll concentrate less on whether that explanation or another makes more sense and concentrate on the empirics–all hail data science!

Data

For this endeavour I collected data from a number of sources–some on wikipedia, others in reports.

Methodology

Constructing measurement models is a hobby of mine.

Pre-processesing

First, let’s scale the data. Scaling is important because variables measured on very different scales will lead to spurious results, due to the larger weight exerted by variables captured in numerically larger values. There are several scaling methods to choose from, but I’ll go with a standard conversion to a standard normal variable (e.g. distributed mean 0, standard deviation 1), through subtracting the mean and dividing by the standard deviation. Note here that I coded variables in such a manner as to presumably achieve positive correlations as opposed to a mix of positive and negative correlations–so sometimes ranges are [0, -inf) as opposed to [0, +inf). This just makes a few steps easier down the road.

Note that you’ll see a subsetting of the data: I initially collected data for more regions than I ultimately opted to include. You’ll also see an imputation step due to limited availability for some measurements. The countries with the largest missingness problems (e.g. Pacific Islands) are excluded from eventual analysis. This leaves 64 countries with 53 total

## 
##  iter imp variable
##   1   1  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   1   2  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   1   3  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   1   4  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   1   5  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   2   1  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   2   2  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   2   3  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   2   4  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   2   5  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   3   1  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   3   2  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   3   3  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   3   4  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   3   5  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   4   1  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   4   2  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   4   3  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   4   4  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   4   5  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   5   1  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   5   2  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   5   3  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   5   4  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI
##   5   5  Religion.Important  PubOp  Church  Suic  Water  COL  Internet  Safety  Specialists  DocSal  DocRelSal  Irreligion  Informal.Payments  Power.Outages  Prosperity.Index  Bribes  Mercer  Innovation  RuleOfLaw  Satisfaction  CapRent  CountryAptCost  SocialCap  AdulteryUnac  HomUn  AbortionOp  PremSexUnac  GAI

Pre-processessing visualization

To get an early hint of the underlying dimensionality of the data, let’s perform a singular value decomposition and visualize the singular values. In general this will resemble an elbow. If there was one latent variable reasonably capturing all of these different measurements, then there would be one large value followed by k-1 small ones, where k is the total number of variables in the data.

svd(df[,c(3:k)])$d %>% plot(pch=16, main="Singular Value Decomposition of Modernization Variables", xlab="Singular Value", ylab="Magnitude of Singular Value")

As another way of visualizing, let’s do a correlation plot.

corrplot(cor(df[,c(3:k)]))

Scale Creation

Ok, now for some fun stuff. And by fun I mean “math that R will do for me.”

Development

When I first envisioned this endeavor, I preferred creating separate scales for governmental development and economic development. However, I quickly realized that separating the two did not make sense. The development of governmental institutions (reflected in metrics such as Rule of Law and Judicial Independence from the Index of Economic Freedom, and World Bank data such as informal payments being made) correlate so strongly with measurements of economic/social development (life expectancy, GDP Per Capita) that concentrating on the two separately would introduce redundancy. It would be redundantly redundant (wink wink).

As such, I examined a reliability measurement (Cronbach’s Alpha) of the selected modernization variables. And as you’ll see below, the alpha is pretty good, well above the informal threshold for what constitutes a reasonable scale. This isn’t exactly the most surprising thing though–is it a shocker that countries with more bribes have less governmental integrity? Some of these are essentially different organizations’ measurements of the same thing. Luckily, while correlated, they diverge sufficiently so as to not introduce mathematical redundancy, which would prevent matrices from being able to be inverted along with other such horrors.

## 
## Reliability analysis   
## Call: psych::alpha(x = df %>% select(Life.Expectancy, GDP.Per.Capita, 
##     Internet, Prosperity.Index, Innovation, RuleOfLaw, Integrity, 
##     Judicial, Pollution, Power.Outages, AirPollDeaths, Informal.Payments, 
##     Bribes, Fragility, SocialCap, HospitalBeds))
## 
##   raw_alpha std.alpha G6(smc) average_r S/N    ase    mean  sd median_r
##       0.96      0.96    0.98      0.62  26 0.0065 3.2e-17 0.8     0.66
## 
##  lower alpha upper     95% confidence boundaries
## 0.95 0.96 0.98 
## 
##  Reliability if an item is dropped:
##                   raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## Life.Expectancy        0.96      0.96    0.98      0.62  24   0.0070 0.055
## GDP.Per.Capita         0.96      0.96    0.98      0.61  23   0.0073 0.055
## Internet               0.96      0.96    0.98      0.63  26   0.0066 0.057
## Prosperity.Index       0.96      0.96    0.98      0.60  22   0.0075 0.051
## Innovation             0.96      0.96    0.98      0.61  23   0.0073 0.054
## RuleOfLaw              0.96      0.96    0.98      0.60  23   0.0075 0.052
## Integrity              0.96      0.96    0.98      0.60  23   0.0075 0.052
## Judicial               0.96      0.96    0.98      0.61  23   0.0072 0.054
## Pollution              0.96      0.96    0.98      0.63  25   0.0068 0.055
## Power.Outages          0.96      0.96    0.98      0.64  27   0.0063 0.056
## AirPollDeaths          0.96      0.96    0.98      0.63  25   0.0067 0.053
## Informal.Payments      0.96      0.96    0.98      0.64  27   0.0064 0.052
## Bribes                 0.96      0.96    0.98      0.61  24   0.0072 0.056
## Fragility              0.96      0.96    0.98      0.60  23   0.0074 0.053
## SocialCap              0.96      0.96    0.98      0.61  24   0.0071 0.057
## HospitalBeds           0.97      0.97    0.98      0.68  32   0.0055 0.025
##                   med.r
## Life.Expectancy    0.65
## GDP.Per.Capita     0.65
## Internet           0.69
## Prosperity.Index   0.65
## Innovation         0.65
## RuleOfLaw          0.65
## Integrity          0.65
## Judicial           0.65
## Pollution          0.66
## Power.Outages      0.69
## AirPollDeaths      0.67
## Informal.Payments  0.69
## Bribes             0.66
## Fragility          0.65
## SocialCap          0.66
## HospitalBeds       0.69
## 
##  Item statistics 
##                    n raw.r std.r r.cor r.drop     mean sd
## Life.Expectancy   64  0.84  0.84  0.84   0.81 -5.8e-17  1
## GDP.Per.Capita    64  0.90  0.90  0.90   0.88  5.7e-19  1
## Internet          64  0.71  0.71  0.69   0.67 -1.4e-16  1
## Prosperity.Index  64  0.97  0.97  0.98   0.97  6.2e-16  1
## Innovation        64  0.91  0.91  0.91   0.89  2.8e-16  1
## RuleOfLaw         64  0.95  0.95  0.96   0.94  6.3e-17  1
## Integrity         64  0.95  0.95  0.96   0.94  1.3e-16  1
## Judicial          64  0.89  0.89  0.89   0.88 -1.4e-16  1
## Pollution         64  0.76  0.76  0.75   0.73 -1.4e-16  1
## Power.Outages     64  0.64  0.64  0.62   0.59  9.0e-18  1
## AirPollDeaths     64  0.73  0.73  0.73   0.69 -8.0e-17  1
## Informal.Payments 64  0.62  0.62  0.59   0.57  8.2e-17  1
## Bribes            64  0.87  0.87  0.86   0.85 -1.6e-17  1
## Fragility         64  0.94  0.94  0.94   0.93  1.9e-17  1
## SocialCap         64  0.86  0.86  0.85   0.84 -2.7e-16  1
## HospitalBeds      64  0.29  0.29  0.25   0.21  8.0e-17  1

Ok, convinced of the merits of my lovely scale? I sure hope so. Onwards! Let’s actually convert this to a single measurement. We’ll use what is probably the most common method around–exploratory factor analysis.

df5 <- df %>% select(
  Life.Expectancy, GDP.Per.Capita,Internet,
  Prosperity.Index,Innovation,RuleOfLaw, Integrity, Judicial,
  Power.Outages,AirPollDeaths,Bribes, Fragility,
  SocialCap, Mercer, Nurses, HospitalBeds
)
row.names(df5) <- df$Country
Development <- factanal(df5, factors=1, scores='regression')
Development
## 
## Call:
## factanal(x = df5, factors = 1, scores = "regression")
## 
## Uniquenesses:
##  Life.Expectancy   GDP.Per.Capita         Internet Prosperity.Index 
##            0.338            0.181            0.536            0.021 
##       Innovation        RuleOfLaw        Integrity         Judicial 
##            0.154            0.059            0.069            0.196 
##    Power.Outages    AirPollDeaths           Bribes        Fragility 
##            0.696            0.504            0.253            0.092 
##        SocialCap           Mercer           Nurses     HospitalBeds 
##            0.293            0.139            0.361            0.948 
## 
## Loadings:
##                  Factor1
## Life.Expectancy  0.814  
## GDP.Per.Capita   0.905  
## Internet         0.681  
## Prosperity.Index 0.989  
## Innovation       0.920  
## RuleOfLaw        0.970  
## Integrity        0.965  
## Judicial         0.897  
## Power.Outages    0.552  
## AirPollDeaths    0.704  
## Bribes           0.864  
## Fragility        0.953  
## SocialCap        0.841  
## Mercer           0.928  
## Nurses           0.800  
## HospitalBeds     0.228  
## 
##                Factor1
## SS loadings     11.158
## Proportion Var   0.697
## 
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 305.96 on 104 degrees of freedom.
## The p-value is 9.15e-22
Development$scores
##                            Factor1
## Albania                -0.88064231
## Argentina              -0.42545636
## Armenia                -0.67876639
## Australia               1.22086443
## Austria                 1.18512426
## Belarus                -0.77869809
## Belgium                 1.03196382
## Belize                 -1.04932109
## Bolivia                -1.47655323
## Bosnia and Herzegovina -0.86088639
## Brazil                 -0.69985438
## Bulgaria               -0.35213578
## Canada                  1.24703984
## Chile                   0.25878189
## Colombia               -0.85041962
## Costa Rica              0.06906005
## Croatia                -0.09581506
## Cyprus                  0.21747533
## Czech Republic          0.53775705
## Denmark                 1.62653447
## Ecuador                -1.00036184
## El Salvador            -1.12527114
## Estonia                 0.91783286
## Finland                 1.30588674
## France                  0.94141051
## Georgia                -0.50080105
## Germany                 1.31351066
## Greece                 -0.05057861
## Guatemala              -1.36071022
## Honduras               -1.47863551
## Hungary                -0.16110829
## Iceland                 1.34551979
## Ireland                 1.18637562
## Italy                   0.42929085
## Latvia                  0.14359776
## Lithuania               0.32928796
## Mexico                 -0.83548300
## Moldova                -1.06405667
## Montenegro             -0.37329942
## Netherlands             1.42089037
## New Zealand             1.38094200
## Nicaragua              -1.51419947
## North Macedonia        -0.72235418
## Norway                  1.64040313
## Panama                 -0.45815580
## Paraguay               -1.02238349
## Peru                   -0.81031461
## Poland                  0.21142643
## Portugal                0.67466322
## Romania                -0.11687965
## Russia                 -0.88867547
## Serbia                 -0.56701237
## Slovakia                0.13489740
## Slovenia                0.63539963
## Spain                   0.74323570
## Suriname               -0.99024716
## Sweden                  1.52899112
## Switzerland             1.57263984
## Ukraine                -0.99415617
## United Kingdom          1.13670182
## United States           0.95362966
## Uruguay                 0.22013860
## Venezuela              -2.22790680
## Guyana                 -1.15013315

Secularism

I next create a scale for secularism. Well, I call it faith, but that’s just because it comes back decreasing in secularism and multiplying by -1 is arbitrarily out of the question. But basically, more secularization = less faith and vice versa. Below I select variables, calculate Cronbach’s Alpha to assess scale reliability, and view the results of the singular value decomposition of the matrix of values as an additional robustness check, because we all need a little more SVD in our lives.

psych::alpha(df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
  AdulteryUnac,HomUn,AbortionOp,PremSexUnac, GAI))
## Number of categories should be increased  in order to count frequencies.
## 
## Reliability analysis   
## Call: psych::alpha(x = df %>% select(Religion.Important, PubOp, Church, 
##     Abor, Irreligion, AdulteryUnac, HomUn, AbortionOp, PremSexUnac, 
##     GAI))
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase     mean   sd median_r
##       0.91      0.91    0.94      0.49 9.8 0.018 -2.7e-17 0.74     0.51
## 
##  lower alpha upper     95% confidence boundaries
## 0.87 0.91 0.94 
## 
##  Reliability if an item is dropped:
##                    raw_alpha std.alpha G6(smc) average_r  S/N alpha se var.r
## Religion.Important      0.89      0.89    0.93      0.48  8.4    0.020 0.035
## PubOp                   0.90      0.90    0.93      0.49  8.7    0.019 0.033
## Church                  0.90      0.90    0.93      0.49  8.6    0.020 0.038
## Abor                    0.91      0.91    0.93      0.52  9.6    0.018 0.034
## Irreligion              0.89      0.89    0.92      0.48  8.1    0.021 0.033
## AdulteryUnac            0.91      0.91    0.95      0.54 10.7    0.016 0.027
## HomUn                   0.89      0.89    0.93      0.47  8.1    0.021 0.036
## AbortionOp              0.89      0.89    0.93      0.47  7.8    0.022 0.034
## PremSexUnac             0.90      0.90    0.94      0.49  8.8    0.020 0.040
## GAI                     0.90      0.90    0.93      0.51  9.2    0.018 0.025
##                    med.r
## Religion.Important  0.51
## PubOp               0.51
## Church              0.48
## Abor                0.55
## Irreligion          0.46
## AdulteryUnac        0.57
## HomUn               0.46
## AbortionOp          0.45
## PremSexUnac         0.51
## GAI                 0.51
## 
##  Item statistics 
##                     n raw.r std.r r.cor r.drop     mean sd
## Religion.Important 64  0.79  0.79  0.79   0.73 -6.5e-19  1
## PubOp              64  0.75  0.75  0.74   0.68  1.2e-18  1
## Church             64  0.76  0.76  0.74   0.70 -8.8e-17  1
## Abor               64  0.63  0.63  0.60   0.54  1.9e-17  1
## Irreligion         64  0.83  0.83  0.83   0.78  5.0e-18  1
## AdulteryUnac       64  0.49  0.49  0.41   0.38 -1.3e-17  1
## HomUn              64  0.83  0.83  0.81   0.78 -9.8e-18  1
## AbortionOp         64  0.87  0.87  0.87   0.84  2.6e-17  1
## PremSexUnac        64  0.74  0.74  0.70   0.67 -6.9e-18  1
## GAI                64  0.68  0.68  0.67   0.59 -2.1e-16  1
svd(df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
                  AdulteryUnac,HomUn,AbortionOp,PremSexUnac))$d %>% 
  plot(pch=16, main="Results of Singular Value Decomposition of Secularization Variables", xlab="Singular Value", ylab="Magnitude of Singular Value")

Ok, Cronbach’s Alpha looks ok (a bit less impressive than the first one but hey it has a lot of human questions in it and humans aren’t necessarily so consistent in their answers), and demonstrated with the SVD that one latent variable captures a good bit of the action.

With that being done, let’s rinse and repeat with some more factor analysis.

Faith <-df %>% select(Religion.Important,PubOp,Church,Abor,Irreligion,
                      AdulteryUnac,HomUn,AbortionOp,PremSexUnac)
row.names(Faith) <- df$Country
Faith <- factanal(Faith, factors=1, scores='regression')
Faith
## 
## Call:
## factanal(x = Faith, factors = 1, scores = "regression")
## 
## Uniquenesses:
## Religion.Important              PubOp             Church               Abor 
##              0.300              0.561              0.511              0.685 
##         Irreligion       AdulteryUnac              HomUn         AbortionOp 
##              0.244              0.857              0.426              0.210 
##        PremSexUnac 
##              0.524 
## 
## Loadings:
##                    Factor1
## Religion.Important 0.837  
## PubOp              0.663  
## Church             0.700  
## Abor               0.561  
## Irreligion         0.870  
## AdulteryUnac       0.378  
## HomUn              0.758  
## AbortionOp         0.889  
## PremSexUnac        0.690  
## 
##                Factor1
## SS loadings      4.683
## Proportion Var   0.520
## 
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 86.75 on 27 degrees of freedom.
## The p-value is 3.45e-08
Faith$scores
##                            Factor1
## Albania                 0.27460352
## Argentina              -0.14545910
## Armenia                 1.39605320
## Australia              -1.05719038
## Austria                -0.85506924
## Belarus                -0.15603821
## Belgium                -1.40247968
## Belize                  0.61591818
## Bolivia                 1.13317521
## Bosnia and Herzegovina  0.69938286
## Brazil                  0.74817131
## Bulgaria               -0.07140885
## Canada                 -1.03093144
## Chile                  -0.04320130
## Colombia                0.87728276
## Costa Rica              0.76598928
## Croatia                 0.42238665
## Cyprus                  0.36849602
## Czech Republic         -1.39661687
## Denmark                -1.53893430
## Ecuador                 1.12393797
## El Salvador             1.52338039
## Estonia                -1.29150759
## Finland                -1.12811617
## France                 -1.47303927
## Georgia                 1.05818497
## Germany                -1.28874922
## Greece                  0.17678240
## Guatemala               1.81152952
## Honduras                0.04048214
## Hungary                -0.93199542
## Iceland                -1.09271411
## Ireland                -1.07336030
## Italy                  -0.13686091
## Latvia                 -0.60043094
## Lithuania              -0.08941334
## Mexico                  0.33703252
## Moldova                 0.77741428
## Montenegro              0.52593919
## Netherlands            -1.20332816
## New Zealand            -0.89076086
## Nicaragua               1.28617523
## North Macedonia         1.13095654
## Norway                 -1.42938737
## Panama                  1.19877090
## Paraguay                1.44781723
## Peru                    0.88068809
## Poland                  0.51477570
## Portugal               -0.10746666
## Romania                 0.65611123
## Russia                 -0.13570222
## Serbia                  0.20964138
## Slovakia                0.35566716
## Slovenia               -0.80563264
## Spain                  -1.12101382
## Suriname                1.60654850
## Sweden                 -1.72328246
## Switzerland            -0.75727603
## Ukraine                 0.41079574
## United Kingdom         -1.27618532
## United States          -0.03895117
## Uruguay                 0.25395288
## Venezuela               0.78406461
## Guyana                  0.88039578

Modernization and Secularization

Having created two scales with sufficient reliability, let’s go back to our questions:

Below, I generate two plots. The first one pools together all countries for a scatterplot with an imposed regression line constituting the best linear fit for faith against development. We see a negative and discernible slope–in other words, that slope surely can’t be argued to be 0. I additionally provide the location of the coordinates of the United States, and a -45 degree line associated with it, from which we can identify which countries are at a comparable level of associated modernization and secularism. The answer to this question is…basically none.

The second plot creates separate regression lines for the different continents (ignoring that I’ve fancifully created a “continent” of U.S./Canada/Australia/New Zealand, all settler societies). We see that all the slopes are negative, though Latin America’s looks a bit flatter than the other ones. In other words, wealthier countries are basically more secular in general, but in Latin America, the relationship seems less strong than in other continents.

A few caveats are in order, of course:

scores <- data.frame(
  Country=df$Country %>% as.vector,
  Continent=df$Continent %>% as.vector,
  Development=Development$scores[,1] %>% as.vector,
  Faith=Faith$scores %>% as.vector
)
xintercept=scores %>% filter(Country=="United States") %>% pull(Development)
yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)
p <- ggplot(scores,
            aes(x=Development, y=Faith,
                label=Country)) + 
  geom_point(aes(colour=Continent)) +
  geom_text_repel() +
  theme_dark() +
  geom_smooth(method='lm') +
  geom_vline(xintercept=scores %>% filter(Country=="United States") %>% pull(Development)) +
  geom_hline(yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)) +
  geom_abline(intercept=xintercept+yintercept,
              slope=-1)
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues
p <- ggplot(scores,
            aes(x=Development, y=Faith,
                label=Country, colour=Continent,
                group=Continent)) + 
  geom_point() +
  geom_text_repel() +
  theme_dark() +
  geom_smooth(method='lm') +
  geom_vline(xintercept=scores %>% filter(Country=="United States") %>% pull(Development)) +
  geom_hline(yintercept=scores %>% filter(Country=="United States") %>% pull(Faith)) +
  geom_abline(intercept=xintercept+yintercept,
              slope=-1)
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'
## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

## Warning in geom2trace.default(dots[[1L]][[4L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

But wait, there’s more…

A simple statistical model

What would such a report be without a statistical model? Well, shorter, for one thing.

We’ll start with a simple model–literally, simple linear regression. To be a bit more proper, we’ll assess what effect a one “unit” change in development has on Faith. Of course, these units are constructs, not empirical quantities, so we might want to change our parlance a bit–does modernization decrease faith by a little, a lot, or not at all? Or might we see something completely unexpected? (Spoilers: No)

mod <- lm(Faith~Development, data=scores)
summary(mod)
## 
## Call:
## lm(formula = Faith ~ Development, data = scores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17245 -0.35131 -0.00739  0.43739  0.83926 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.912e-17  6.491e-02    0.00        1    
## Development -8.203e-01  6.570e-02  -12.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5193 on 62 degrees of freedom
## Multiple R-squared:  0.7155, Adjusted R-squared:  0.7109 
## F-statistic: 155.9 on 1 and 62 DF,  p-value: < 2.2e-16
plot_model(mod, type='eff')
## $Development

Survey says….surprise surprise, modernization has a negative effect on faith. A shocking, breathtaking conclusion says…nobody. The slope is certainly something to type home about. To put it in perspective, let’s examine our dependent variable:

summary(scores$Faith)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.7233 -0.9567  0.1086  0.0000  0.7688  1.8115

The minimum value is -1.7232825, and the maximum value is 1.8115295. Some fancy math tells us the range of the variable is -3.534812. The coefficient for the model tells us that going one unit up in Modernization decreases that Faith variable by -0.8203049 (thank you, R weirdness with single vs. double brackets). So, think about how much relatively you go down in faith with a single unit change in Modernization.

Ok, I’ve got a fever, and the only cure is Tylenol. No, actually, it’s…more statistical models. We saw a difference in regions in the plots above so let’s add that to the model: maybe each region has its own baseline.

A slightly more complicated model

mod2 <- lm(Faith~Development + Continent, data=scores)
summary(mod2)
## 
## Call:
## lm(formula = Faith ~ Development + Continent, data = scores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.18519 -0.29383  0.03051  0.45954  0.97720 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        -0.01171    0.10863  -0.108   0.9145    
## Development                        -0.63432    0.12590  -5.038 4.74e-06 ***
## ContinentLatin America              0.29945    0.17455   1.716   0.0915 .  
## ContinentNorth America/Australasia  0.01883    0.33021   0.057   0.9547    
## ContinentWestern Europe            -0.33219    0.24488  -1.357   0.1801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5097 on 59 degrees of freedom
## Multiple R-squared:  0.7391, Adjusted R-squared:  0.7214 
## F-statistic: 41.79 on 4 and 59 DF,  p-value: < 2.2e-16
plot_model(mod2, type='eff')
## $Development

## 
## $Continent

As our plots above had suggested, Latin America seems to be a bit higher in faith than the other three regions, which are otherwise not able to be distinguished one from another. The slope of our Development variable does change a bit–some of its effect in the simple model likely owed to Latin American countries scoring lower on Development while also having lower levels of Secularization.

Are we sure e like the more complicated model? Well, beyond the R squared difference, ANOVA can do a bit for us:

anova(mod2, mod)
## Analysis of Variance Table
## 
## Model 1: Faith ~ Development + Continent
## Model 2: Faith ~ Development
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     59 15.326                           
## 2     62 16.717 -3   -1.3911 1.7851 0.1598

One last little quibble–we may have heteroscedasticity, one of the deadly sins of ordinary least squares estimation. Just in case, since our Breusch-Pagan test is basically borderline, we’ll use some robust standard errors to make ourselves a little happier. I’m happier already, aren’t you?

bptest(mod2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 6.3268, df = 4, p-value = 0.176
coeftest(mod2, vcov = sandwich::vcovHC, type = "HC0")
## 
## t test of coefficients:
## 
##                                     Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)                        -0.011707   0.124850 -0.0938   0.92561    
## Development                        -0.634321   0.141371 -4.4869 3.402e-05 ***
## ContinentLatin America              0.299452   0.161438  1.8549   0.06861 .  
## ContinentNorth America/Australasia  0.018826   0.305046  0.0617   0.95100    
## ContinentWestern Europe            -0.332195   0.268071 -1.2392   0.22018    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A more than slightly more complicated model

The results are in and we are keeping the more complicated model.

But before we go let’s try one more thing–what we did above is called a random intercept model, if we want to get a bit particular. What if we want to test random slope and random intercept? Well, we can do so fairly easily–we’ll add an interaction term for it.

mod3 <- lm(Faith~Development*Continent, data=scores)
summary(mod3)
## 
## Call:
## lm(formula = Faith ~ Development * Continent, data = scores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03934 -0.34304  0.03864  0.29434  0.90003 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                    -0.056211   0.113302  -0.496
## Development                                    -0.813588   0.188429  -4.318
## ContinentLatin America                          0.568494   0.232673   2.443
## ContinentNorth America/Australasia              2.054041   1.977623   1.039
## ContinentWestern Europe                        -0.064857   0.469826  -0.138
## Development:ContinentLatin America              0.429762   0.266461   1.613
## Development:ContinentNorth America/Australasia -1.478803   1.641780  -0.901
## Development:ContinentWestern Europe            -0.007565   0.412825  -0.018
##                                                Pr(>|t|)    
## (Intercept)                                      0.6218    
## Development                                    6.49e-05 ***
## ContinentLatin America                           0.0177 *  
## ContinentNorth America/Australasia               0.3034    
## ContinentWestern Europe                          0.8907    
## Development:ContinentLatin America               0.1124    
## Development:ContinentNorth America/Australasia   0.3716    
## Development:ContinentWestern Europe              0.9854    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5055 on 56 degrees of freedom
## Multiple R-squared:  0.7564, Adjusted R-squared:  0.7259 
## F-statistic: 24.84 on 7 and 56 DF,  p-value: 4.998e-15
plot_model(mod3)

The results are in! And, I don’t think I’m too impressed. Maybe that slope for North America/Oceania is steeper than the rest. But…maybe not? Meh. I think relative simplicity rules on this one.

Conclusions

Let’s turn to the bullet points.

Open Questions